#Пакеты

library(tidyverse)
library(readxl)
library(openxlsx)
library(dplyr)
library(ggplot2)
library(plotly)
library(ggpubr)
library(nortest)
library(skimr)
library(flextable)
library(ggbeeswarm)
library(rstatix)
library(corrplot)
library(corrr)
library(GGally)
library(factoextra)
library(pheatmap)
library(ggbiplot)
#library(tidymodels)
library(embed)
library(FactoMineR)

#Задание 1. Загрузите датасет life_expectancy_data.RDS

readRDS("life_expectancy_data.RDS") -> data_raw
data_raw %>%
  glimpse()

#Задание 2. Сделайте интерактивный plotly график любых двух нумерических колонок. Раскрасть по колонке континента, на котором расположена страна

plot_ly(
  data = data_raw[(data_raw$`Life expectancy`!= 0) & (data_raw$`Tuberculosis Incidence` != 0),],
  x = ~ `Life expectancy`,
  y = ~ `Tuberculosis Incidence`,
  color = ~continent,
  marker = list(
    size = 9,
    line = list(color = 'rgba(152, 0, 0, .6)', 
                width = 1)
  )
)   %>%
  layout(
    title = 'Отношение Ожидаемой продолжительности жизни и Заболеваемости туберкулезом',
    yaxis = list(title = 'Ожидаемая продолжительность жизни',
                 zeroline = FALSE),  
    xaxis = list(title = 'Заболеваемость туберкулезом',
                 zeroline = FALSE))

Задание 3. Проведите тест, на сравнение распределений колонки Life expectancy между группами стран Африки и Америки. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку rstatix.

data_raw %>%
  select('Life expectancy', 'continent') %>%
  filter(continent == "Africa" | continent == "Americas") -> data_life

t_test(data = data_life, formula = `Life expectancy` ~ continent)
## # A tibble: 1 × 8
##   .y.             group1 group2      n1    n2 statistic    df        p
## * <chr>           <chr>  <chr>    <int> <int>     <dbl> <dbl>    <dbl>
## 1 Life expectancy Africa Americas    52    38     -12.3  85.8 1.31e-20
stat.test <- data_life %>%
  group_by(continent) %>%
  t_test(data = data_life, formula = `Life expectancy` ~ continent) %>%
  adjust_pvalue(method = "bonferroni") %>%
  add_significance("p.adj")
stat.test
## # A tibble: 1 × 10
##   .y.   group1 group2    n1    n2 statistic    df        p    p.adj p.adj.signif
##   <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl>    <dbl>    <dbl> <chr>       
## 1 Life… Africa Ameri…    52    38     -12.3  85.8 1.31e-20 1.31e-20 ****
bxp <- ggboxplot(
  data_life, y = "Life expectancy", x = "continent", 
  palette = c("#00AFBB", "#E7B800")
  )

stat.test <- stat.test %>%
  add_xy_position(x = "Life expectancy", dodge = 0.8)
bxp + stat_pvalue_manual(
  stat.test,  label = "p", tip.length = 0
  )

Задание 4. Сделайте новый датафрейм, в котором оставите все численные колонки кроме Year. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

data_selected <- data_raw %>% select(!Year) %>% select(where(is.numeric))

data_cor <- cor(data_selected)

#график 1
corrplot(data_cor, method = 'color', type = "lower")

#график 2
data_cor %>% 
  rplot()

Задание 5. Постройте иерархическую кластеризацию на этом датафрейме.

data_selected_scaled <- scale(data_selected)

data_selected_dist <- dist(data_selected_scaled, 
                        method = "euclidean"
                        )
as.matrix(data_selected_dist)[1:6,1:6]
##          1        2        3        4        5        6
## 1 0.000000 7.605708 6.331840 4.414874 6.645623 7.923487
## 2 7.605708 0.000000 2.624659 7.921597 3.357361 3.631018
## 3 6.331840 2.624659 0.000000 6.321666 4.350331 3.464837
## 4 4.414874 7.921597 6.321666 0.000000 8.095849 7.161240
## 5 6.645623 3.357361 4.350331 8.095849 0.000000 4.966244
## 6 7.923487 3.631018 3.464837 7.161240 4.966244 0.000000
data_selected_hc <- hclust(d = data_selected_dist, 
                        method = "ward.D2")

fviz_dend(data_selected_hc, 
          cex = 0.1)

# Задание 6. Сделайте одновременный график heatmap и иерархической кластеризации. Содержательно интерпретируйте результат.

pheatmap(data_selected_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = data_selected_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5,
         cutree_cols = length(colnames(data_selected_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

#Интерпретация: 
# В первой группе пациентов выделяются показатели смертности и заболеваемости туберкулезом и в отрицательную сторону вакцинация. В целом данные и особенно Пятая группа довольно размыты. Выделилась третья маленькая группа, в ней сильно выделились показатели GDP и GNI.

Задание 7. Проведите PCA анализ на этих данных. Проинтерпретируйте результат.

data_full.pca <- prcomp(data_selected, 
                        scale = T)
summary(data_full.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18      PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 7.273e-16
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.000e+00
# Смотрим Cumulative Proportion - PC1, PC2 и PC3 чуть больше 60% (заметно меньше 75%), для первых двух - 51: (меньше 60%). В целом, это не слишком хороший результат. Далее посомтрим график
fviz_eig(data_full.pca, addlabels = T, ylim = c(0, 40))

# Выделяется первая компонента, вторая и третья заметно меньше. Дальше имеет смысл рассматривать в первую очередь их
# График для двух главных компонент
fviz_pca_var(data_full.pca, col.var = "contrib")

# График для двух главных компонент с выбранным количеством контрибьюторов
fviz_pca_var(data_full.pca, 
             select.var = list(contrib = 8), 
             col.var = "contrib")

# Стрелки довольно близки к кругу
fviz_contrib(data_full.pca, choice = "var", axes = 1, top = 24) # 1

fviz_contrib(data_full.pca, choice = "var", axes = 2, top = 24) # 2

fviz_contrib(data_full.pca, choice = "var", axes = 3, top = 24) # 3

# Из чего состоят три главные компоненты
#PC1 - довольно много переменных, основных - 4
#PC2 составляют в первую очередь переменные иммунизаций
#PC3 состоит из GDP и GNI

Задание 8. Постройте biplot график для PCA. Раскрасьте его по значениям континентов. Переведите его в plotly. Желательно, чтобы при наведении на точку, вы могли видеть название страны.

Вариант 1

ggbiplot(data_full.pca, 
         scale=0, alpha = 0.1, varname.size = 4, groups = data_raw$continent, labels = data_raw$Country, labels.size = 1, ellipse = T) + 
  theme_minimal() -> PCA_biplot

ggplotly(PCA_biplot, tooltip = c("groups", "labels"))

Вариант 2

ggbiplot(data_full.pca, 
         scale=0, alpha = 0.3, groups = data_raw$continent, ellipse = T) + 
  theme_minimal() -> PCA_biplot

ggplotly(PCA_biplot, tooltip = c("groups"))

Задание 9. Дайте содержательную интерпретацию PCA анализу.

Выделяются 3 главные компоненты. В первую компоненту входят сразу несколько разных показателей (основная - ожидаемая продолжительность жизни), во вторую - показатели иммунизации (они хорошо скоррелированы), третья компонента состоит из GDP и GNI и они сильно скоррелированы, хотя их вклад заметно меньше других переменных (кроме суицидов, безработицы и лечения туберкулеза) 

Такие показатели, как Городское население и Сельское население, противоположно направлены, но вносят одинаковый вклад. При этом показатель Per Capita совпадает с Urban population по направлению, но вносит несколько меньший вклад. 

Для данного датасета имеет смысл проводить PCA, но

Задание 10. Сравните результаты отображения точек между алгоритмами PCA и UMAP.

umap_prep <- recipe(~., data = data_selected) %>% 
  step_normalize(all_predictors()) %>% 
  step_umap(all_predictors()) %>% 
  prep() %>% 
  juice() 
umap_prep2 <- cbind(umap_prep, data_raw)

umap_prep2 %>%
  ggplot(aes(UMAP1, UMAP2)) + 
  geom_point(aes(color = continent,
             alpha = 0.6, size = 2)) +
  labs(color = NULL) 

В UMAP прослеживается тенденция, что данные внутри региона близки друг к другу, тк на получившемся графике точки группируются по регионам. Для регионов Азия, Америка и Океания отмечатеся более разраженное расположение точек, что возможно объясняется большец неоднорожностью данных внутри региона.